Computing power series expansions of modular 
forms 
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Abstract We exhibit a method to numerically compute power series expansions 
of modular forms on a cocompact Fuchsian group, using the explicit computation a 
fundamental domain and linear algebra. As applications, we compute Shimura curve 
parametrizations of elliptic curves over a totally real field, including the image of 
CM points, and equations for Shimura curves. 



1 Introduction 



A classical modular form / : "K — > C on the upper half -plane "K satisfies the trans- 
lation invariance f(z + 1) = f(z) for z <G "K, so / admits a Fourier expansion (or 
q-expansion) 

oo 

n=0 

at the cusp oo, where q = e 2mz . If further / is a normalized eigenform for the Hecke 
operators T n , then the coefficients a„ are the eigenvalues of T„ for n relatively prime 
to the level of /. The q-expansion principle expresses in a rigorous way the fact that 
a modular form is characterized by its ^-expansion, and for this reason (and others) 
^-expansions remain an invaluable tool in the study of classical modular forms. 

By contrast, modular forms on cocompact Fuchsian groups do not admit q- 
expansions due to the lack of cusps. A modular form f : — > C of weight k e 2Z>o 
for a cocompact Fuchsian group F < PSL.2(M) is a holomorphic map satisfying 
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f(gz)=j(g,z) k f(z) 



for all gef, where j{g,z) = cz + d if g ■ 



a b 
c d 



However, not all is lost: such a modular form / still admits a power series ex- 
pansion in the neighborhood of a point p e !K. Indeed, a ^-expansion is really just 
a power series expansion at °° in the parameter q, convergent for \q\ < 1; so it is 
natural to consider a neighborhood of p normalized so the expansion also converges 
in the unit disc but for a parameter w. For this purpose, we map the upper half -plane 
conformally to the unit disc D via the map 

z-p 



z w(z) 



z-p 



sending p ^ w(p) = 0, where denotes complex conjugation. We consider the 
series expansion of a form / of weight k given by 

oo 

/(z) = (l-w)*£>w n (*) 

n=0 

where w = w(z), convergent in the disc D with \w\ < 1. 

There are several reasons to consider series of the form (*). First, the term 

l-w(z) = ^ 
z-p 

is natural to include as it arises from the automorphy factor of the linear fractional 
transformation w(z). Second, the ordinary Taylor coefficients arise from evaluating 
derivatives of /, but the derivative of a modular form of weight k ^ is no longer 
a modular form! This can be ameliorated by considering instead the differential 
operator introduced by Maass and Shimura: 

1 (d k 



2ni \dz z — z 

If / is a modular form of weight k then d^f transforms like a modular form of 
weight k + 2, but at the price that / is now only real analytic. The coefficients b n in 
the expansion (*) then arise from evaluating Shimura-Maass derivatives d"f of / 
at the point p, where we let d n = d k+2 ( n -\) o ... o ° dk- Finally, when p is a CM 
point and r is a congruence group, the coefficients b n are algebraic (up to a nor- 
malization factor), as shown by Shimura [36, 37]. Rodriguez-Villegas and Zagier 
[34] and O' Sullivan and Risager [33] (see also the exposition by Zagier [48] and 
further work by Bertolini, Darmon, and Prasanna [6, §§5-6]) further linked the coef- 
ficients b n to square roots of central values of Rankin-Selberg L-series and show that 
they satisfy a recursive formula arising from the differential structure of the ring of 
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modular forms. Mori [29, 30] and Datskovsky and Guerzhoy [14] have also studied 
the /?-adic interpolation properties of the coefficients b n in the context of p-adic L- 
functions, with some results also in the cocompact case. For these reasons, we will 
consider series expansions of the form (*) in this article, which we will call power 
series expansions. 

In the special case when the Fuchsian group F is commensurable with a triangle 
group, the differential approach to power series expansions of modular forms is 
particularly pleasing. For a,b,c e Z>2 U {°°} with l/a+l/b + l/c < l,we define 
the (fl,fo,c)-triangle group to be the subgroup of orientation-preserving isometries 
in the group generated by reflections in the sides of a hyperbolic triangle with angles 
7t/a,n/b,n/c. In this case, power series expansions for a uniformizing function at 
the vertices of the fundamental triangle are obtained as the inverse of the ratio of 
2-Fi-hypergeometric functions. This case was also of great classical interest, and has 
been taken up again more recently by Bayer [2], Bayer and Travesa [3, 4], the first 
author [45], and Baba and Granath [1]; this includes the well-studied case where 
the Fuchsian group arises from the quaternion algebra of discriminant 6 over Q, 
corresponding to the (2,4,6)-triangle group. 

In this article, we exhibit a general method for numerically computing power 
series expansions of modular forms for cocompact Fuchsian groups. It is inspired by 
the method of Stark [41] and Hejhal [22, 23, 24], who used the same basic principle 
to compute Fourier expansions for Maass forms on SL.2(Z) and the Hecke triangle 
groups. (There has been substantial further work in this area; see for example Then 
[43], Booker, Strombergsson, and Venkatesh [8], and the references therein.) 

The basic idea is quite simple. Let F be a Fuchsian group with compact funda- 
mental domain D C D contained in a circle of radius p > 0. To find a power series 
expansion (*) for a form / on F of weight k valid in D to some precision e > 0, we 
consider the approximation 

f(z)^f N (z) = (l-w) k tbnW" 

n=0 

valid for \w\ < p and some N € Z>o- Then for a point w = w(z) on the circle of 
radius p with w £ D, we find g e F such that w 1 = gw € D; letting z' = z(w r ), by the 
modularity of / we have 

(1 - w'f £ b n (w')" = f N (z') « f(z') = j(g,z) k f(z) « j(g,z) k (l wf £ b n w" 

n=0 n=0 

valid to precision e > 0. For each such point w, this equality imposes an approximate 
(nontrivial) linear relation on the unknown coefficients b„. By taking appropriate 
linear combinations of these relations, or (what seems better in practice) using the 
Cauchy integral formula, we recover the approximate coefficients b„ using standard 
techniques in linear algebra. 

An important issue that we try to address in this paper is the numerical stability 
of this method — unfortunately, we cannot benefit from the exponential decay in the 
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terms of a Fourier expansion as in previous work, so our numerical methods must be 
correspondingly more robust. Although we cannot prove that our results are correct, 
there are several tests that allow one to be quite convinced that they are correct, and 
we show in several examples that they agree with cases that are known. (See also 
Remark 3.2 and work of Booker, Strombergsson and Venkatesh [8], who rigorously 
verify the numerical computations in Hejhal's method for Maass forms.) 

Nelson [31, 32] finds power series expansions by directly computing the Shimizu 
lift (a realization of the Jacquet-Langlands correspondence) of a modular form on 
a Shimura curve over Q to a classical modular curve. It will be interesting to com- 
pare the two techniques. Owing to its simplicity, we believe that our method is 
worthy of investigation. It has generalizations to a wide variety of settings: non- 
congruence groups, nonarithmetic groups (e.g. nonarithmetic triangle groups), real 
analytic modular forms, and higher dimensional groups; and it applies equally well 
for arithmetic Fuchsian groups with an arbitrary totally real trace field. 

Finally, this analysis at the complex place suggest that p-adic analytic methods 
for computing power series expansions would also be interesting to investigate. For 
example, Franc [19] investigates the values of Shimura-Maass derivatives of modu- 
lar forms at CM points from a rigid analytic perspective. 

This paper is organized as follows. In Section 2, we introduce some basic nota- 
tion and background. Then in Section 3, we exhibit our algorithm in detail. Finally, 
in Section 4 we give several examples: in the first two examples we verify the cor- 
rectness of our algorithm; in the third example, we compute the Shimura curve 
parametrization of an elliptic curve over a totally real field and the image of CM 
point; in the fourth example, we show how our methods can be used to compute the 
equation of a Shimura curve. 



2 Preliminaries 

We begin by considering the basic setup of our algorithm; as basic references we 
refer to Beardon [5] and Katok [25]. 



Fuchsian groups and fundamental domains 

Let r be a Fuchsian group, a discrete subgroup of PSL2(M), the orientation- 
preserving isometries of the upper half-plane "K. Suppose that F is cofinite, so 
X = F\ JC has finite hyperbolic area; then F is finitely generated. If X is not com- 
pact, then X has a cusp, and the existence of ^-expansions at cusps in many cases 
obviates the need to compute power series expansions separately (and for Maass 
forms and generalizations, we refer to the method of Hejhal [23]); more seriously, 
our method apparently does not work as well in the non-cocompact case (see Ex- 
ample 2 in Section 4). So we suppose that F is cocompact. 
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We will frequently move between the upper half-plane model K and the Poincare 
unit disc model D for hyperbolic space, which are conformally identified via the 
maps 

w:K -^D z:D -> IK 

/ n z-P i \ pw-p (1) 

z^w(p;z) = — - w ^ z(p;w) = f- 

z—p w—1 

for a choice of point p £ IK. Via this identification, the group F acts also on 2), and 
to ease notation we identify these actions. 

A Fuchsian group F is exacf if there is a finite set G C SI^-ST) with iiT^QHR 
a number field whose image in PSL2(/T) C PSL2(M) generates F. When speaking 
in an algorithmic context, we will suppose that the group F is exact. (Even up to 
conjugation in PSL2(M), not every finitely generated Fuchsian group is exact.) Al- 
gorithms for efficiently computing with algebraic number fields are well-known (see 
e.g. Cohen [12]), and we will use these without further mention. 

In this setting, we have the following result, due to the first author [47]. For a 
point p G IK, we denote by T p = {g £ F : g(p) — p} the stabilizer of p in F. 

Theorem 2.1. There exists an algorithm that, given as input an exact, cocompact 
Fuchsian group F and a point p £ K with F ; , = { 1 }, computes as output a funda- 
mental domain D(p) C K for F and an algorithm that, given z £ K returns a point 
z' £ D(p) and g £ F such that z' — gz. 

In particular, in the course of exhibiting this algorithm a suite of methods for 
computing in the unit disc are developed. The algorithm in Theorem 2. 1 has been 
implemented in Magma [9]. 

The fundamental domain D{p) in Theorem 2.1 is the Dirichlet domain centered 
at p, 

D{p) = {z e K : d(z,p) < d(gz,p) for all g e F} 

where d is the hyperbolic distance. The set D{p) is a closed, connected, and hyper- 
bolically convex domain whose boundary consists of finitely many geodesic seg- 
ments. The image of D(p) in D is analogously described, and we use the same 
notation for it. A domain with this description is indeed desirable for consideration 
of power series centered at p, as we collect in each orbit of F the points closest to 
P- 



Modular forms and power series expansions 

A (holomorphic) modular form of weight k £ 2Z>o for F is a holomorphic map 
/ : K -> C with the property that 



f(gz) = j(g,z) k f(z) 



(2) 



6 



John Voight and John Willis 



for all g = ± e F where j(g,z) = cz + d. (Note that although the matrix is 

only defined up to sign, the expression j(g,z) k is well-defined since k is even. Our 
methods would extend in a natural way to forms of odd weight with character on 
subgroups of SL2QR), but for simplicity we do not consider them here.) Let M^(F) 
be the finite-dimensional C-vector space of modular forms of weight k for F, and 
let M(r) = @kMk(r) be the ring of modular forms for F under multiplication. 
A function / : Jf — > C is said to be nearly holomorphic if it is of the form 



d=Q 

where each : Ji — > C is holomorphic. A nearly holomorphic function is real an- 
alytic. A nearly holomorphic modular form of weight k € 2Z>o for F is a nearly 
holomorphic function / : Ji — > C that transforms under F as in (2). Let A/|(F) be 
the C-vector space of nearly holomorphic modular forms of weight k for F and let 
M*(F) - © fe M|(F); then M*(F) CK* (F). 

Via the identification of the upper half -plane with the unit disc, we can consider 
a modular form / also as a function on the unit disc. The transformation property 
(2) of g e F is then described by 

f{gz) = {cz + dfm = ( {CP+d) JS l {P+d) ) k Az) wherez = z(w). (3) 

Let / be a modular form of weight k for F. Since / is holomorphic in D, it has a 

power series expansion 

f(z) = (l-w) k f J b n w n (4) 

n=0 

with fc„eC and w = w(z), convergent in D and uniformly convergent on any com- 
pact subset. 

The coefficients b n , like the usual Taylor coefficients, are related to derivatives 
of / as follows. Define the Shimura-Maass differential operator : M^(F) -4- 
^ +2 (F)by 

1 (d k 



2ni \dz z — z 

Even if / € M^(F) is holomorphic, d%f € M| +2 (F) is only nearly holomorphic; 
indeed, the ring M* (F) of nearly holomorphic modular forms is the smallest ring of 
functions which contains the ring of holomorphic forms M(F) and is closed under 
the Shimura-Maass differential operators. For n > let dg = dfc+2(n-i) o ■ ■ ■ o d^yi o 
dk, and abbreviate d n = dg. We have the following proposition, proven by induction: 
see e.g. Zagier [48, Proposition 17] (but note the sign error!). 



Lemma 2.2. Let f : JC — > C fee holomorphic at p e M. 77zen / admits the series 
expansion (4) ;'n a neighborhood of p with 
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where y — lia(p). 

Remark 2.3. This expression of the coefficients b„ in terms of derivatives implies 
that they can also be given as (essentially) the constant terms of a sequence of poly- 
nomials satisfying a recurrence relation, arising from the differential structure on 
M*(r): see Rodriguez- Villegas and Zagier [34, §§6-7] and Zagier [48, §5.4, §6.3], 
who carry this out for finite index subgroups of SL2(Z), and more recent work of 
Baba and Granath [1]. 

The expression for the regular derivative 



(2ni) n dz n 

in terms of d n is given in terms of Laguerre polynomials: see Rodriguez- Villegas 
and Zagier [34, §2], Zagier [48, (56)-(57)], or O'Sullivan and Risager [33, Proposi- 
tion 3.1]. 

Lemma 2.4. We have 

"/ n \ (fc + r) n _ r 

and 

D"f = T (-!)""' ( H ) d r f 
r =o \r) (-47cy)"-' 

where y — Imp and (a) m = a(a + 1) • • • (a + m — 1) is the Pochhammer symbol. 



Arithmetic Fuchsian groups 

Among the Fuchsian groups we consider, of particular interest are the arithmetic 
Fuchsian groups. A basic reference is Vigneras [44]; see also work of the first author 
[46] for an algorithmic perspective. 

Let F be a number field with ring of integers 7L F . A quaternion algebra B over F 
is an F-algebra with generators a,j3 € B such that 

a 2 = a, P 2 = b, j3a = ~a[5 

with a,b eF x ; such an algebra is denoted B — 

Let B be a quaternion algebra over F. Then B has a unique (anti-)involution 
~ : B — > B such that the reduced norm nrd(y) = yy belongs to F for all y e B. A 
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place v of F is split or ramified according as B v = B ®p F v = M2(F V ) or not, where 
F v denotes the completion at v. The set S of ramified places of B is finite and of even 
cardinality, and the product D of all finite ramified places is called the discriminant 
ofB. 

Now suppose that F is a totally real field and that B has a unique split real place 
v ^ S corresponding to !«, : B B ® F v = M 2 (R). An orafer C B is a subring 
with FO = B that is finitely generated as a Zf-submodule. Let C B be an or- 
der and let 0| denote the group of units of reduced norm 1 in 0. Then the group 
r B (l) = ioo((V{±l}) C PSL 2 (M) is a Fuchsian group [25, §§5.2-5.3]. An arith- 
metic Fuchsian group T is a Fuchsian group commensurable with r B {\) for some 
choice of B. One can, for instance, recover the usual modular groups in this way, 
taking F = Q, = M 2 (Z) C M 2 (Q) = B, and T C PSL 2 (Z) a subgroup of finite 
index. An arithmetic Fuchsian group T is cofinite and even cocompact, as long as 
B ^ M 2 (Q), which we will further assume. In particular, the fundamental domain 
algorithm of Theorem 2.1 applies. 

Let 9t be an ideal of Zp coprime to D. Define 

0(^)i X = {ye 0:7=1 (mod TO))} 

and let r B (m) = lo<,(0(9T)f ). A Fuchsian group T commensurable with is 
congruence if it contains r B (yX) for some 91. 

The space M^T) has an action of Hecke operators T p indexed by the prime ide- 
als p \ 2)91 that belong to the principal class in the narrow class group of Zp. (More 
generally, one must consider a direct sum of such spaces indexed by the narrow 
class group of F .) The Hecke operators can be understood as averaging over sublat- 
tices, via correspondences, or in terms of double cosets; for a detailed algorithmic 
discussion in this context and further references, see work of Greenberg and the first 
author [21] and Dembele and the first author [16]. The operators T p are semisimple 
and pairwise commute so there exists a basis of simultaneous eigenforms in M^(r) 
for all T p . 

Let K be a totally imaginary quadratic extension of F that embeds in B, and let 
veBbe such that F(v) = K. Let p e 'K be a fixed point of Uo(v). Then we say p 
is a CM point for K. 

Theorem 2.5. There exists Q € C x such that for every CM point pfor K, every con- 
gruence subgroup r commensurable with F B (1), and every eigenform f G M^{T) 
with f(p) £ Q X , we have 

for all n € Z>o- 

The work of Shimura on this general theme of arithmetic aspects of analytic 
functions is vast. For a "raw and elementary" presentation, see his work [36, Main 
Theorem I] on the algebraicity of CM points of derivatives of Hilbert modular forms 
with algebraic Fourier coefficients. He then developed this theme with its connection 
to theta functions by a more detailed investigation of derivatives of theta functions 
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[37]. For the statement above, see his further work on the relationship between the 
periods of abelian varieties with complex multiplications and the derivatives of au- 
tomorphic forms [38, Theorem 7.6]. However, these few references barely scratch 
the surface of this work, and we refer the interested reader to Shimura's collected 
works [40] for a more complete treatment. 

Note that Q. is unique only up to multiplication by an element of Q x . When F = 
Q, the period £2 = £2% can be taken as the Chowla-Selberg period associated to the 
imaginary quadratic field K, given in terms of a product of values of the F -function 
[48, (97)]. (The Chowla-Selberg theory becomes much more complicated over a 
general totally real field; see e.g. Moreno [28].) In the case where ^-expansions are 
available, one typically normalizes the form / to have algebraic Fourier coefficients, 
in which case f(p) € £2 k Q; for modular forms on Shimura curves, it seems natural 
instead to normalize the function so that its values at CM points for K are algebraic. 

For the rest of this section, suppose that F is a congruence arithmetic Fuchsian 
group containing F B (9T). Then from Lemma 2.2, we see that if / e M^(F) and as 
in Theorem 2.5, then the power series (4) at a CM point p can be rewritten as 

f(z)=f(p)(l-w) k Z^(0wy (5) 

«=o 



where 



e = ^mP>* (6) 



fip) V 



and 

= (d n f)(p) ( fip) V = 1 b n (b 
Cn fip) Wf)ip)) n\bo\h 



(7) 



are algebraic. In Section 3, when these hypotheses apply, we will use this extra 
condition to verify that our results are correct. 

Remark 2.6. Returning to Remark 2.3, for a congruence group F, it follows that 
once we have computed just a finite number of the coefficients c n for a finite gen- 
erating set of M(r) to large enough precision to recognize them algebraically and 
thus exactly, following Zagier we can compute the differential structure of M(F) 
and thereby compute exactly all coefficients of all forms in M(F) rapidly. It would 
be very interesting, and quite powerful, to carry out this idea in a general context. 



3 Numerical method 



In this section, we exhibit our method to compute the finite-dimensional space of 
(holomorphic) modular forms M^(F) of weight k e 2Z>o for a cocompact Fuchsian 
group F. Since the only forms of weight are constant, we assume k > 0. 

Remark 3.1. To obtain meromorphic modular forms of weight (for example) we 
can take ratios of holomorphic modular forms of the same weight k. Or, alter- 
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natively, since the Shimura-Maass derivative of a meromorphic modular form of 
weight is in fact meromorphic of weight 2 (k = in the definition), it follows that 
the antiderivative of a meromorphic modular form of weight 2 is meromorphic of 
weight 0. 

Let £ > be the desired precision in the numerical approximations. We write 
x w y to mean approximately equal (to precision e); we leave the task of making 
these approximations rigorous as future work. 

Let p £ "K satisfy r p = {1} and let D = D(p) C D be the Dirichlet domain cen- 
tered at p as in Theorem 2.1, equipped with an algorithm that given a point iveD 
computes geT such that w' — gw € D. Let 

p = p(D) = max{|w| : w G D} 

be the radius of D. 



Approximation by a polynomial 

The series expansion (4) for a form / e Mk(r) converges in the unit disc D and its 
radius of convergence is 1 . To estimate the degree of a polynomial approximation 
for / valid in D to precision e, we need to estimate the sizes of the coefficients b n . 

We observe in some experiments that the coefficients are roughly bounded. If we 
take this as a heuristic, assuming \b n \ < 1 for all n, then the approximation 

Af 

f{z)^f N {z) = {l-wf%b n w n 

n=0 



is valid for all \w\ < p with 



N = 



log(e) 
logp 



(8) 



Remark 3.2. As mentioned in the introduction, if p is a CM point and / is an eigen- 
form for a congruence group, then the coefficient b n is related to the central critical 
value of the Rankin-Selberg L-function L(s,f x 0"), where 9 is the modular form 
associated to a Hecke character for the CM extension given by p. Therefore, the best 
possible bounds on these coefficients will involve (sub)convexity bounds for these 
central L-values (in addition to estimates for the other explicit factors that appear). 

In practice, we simply increase in two successive runs and compare the results 
in order to be convinced of their accuracy. Indeed, a posteriori, we can go back using 
the computed coefficients to give a better estimate on their growth and normalize 
them so that the boundedness assumption is valid. In this way, we only use this 
guess for N in a weak way. 
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The basic idea now is to use the modularity of /, given by (2), at points inside and 
outside the fundamental domain to obtain relations on the coefficients b n . 

For any point iveB with \w\ < p and z = z{w), we can compute gef such that 
w' = gw G D; letting z' = z(w'), by the modularity of / we have 

f N (z') « f(z') = j{g,zff{z) « j(g,z) k f N (z) 

so 

(1 - w'f £ fo„(w')" » 7(«,z)*(l - "0* I b nW " (9) 

n=0 n=0 

and thus 

£*„»&„ «0 (10) 

n=0 

where ("a" for automorphy) 

K H = j(g,z) k (l - w) V - (1 - vw')*(w')". (11) 

If w does not belong to the fundamental domain D, so that g ^ 1, then this imposes 
a nontrivial relation on the coefficients With enough relations, we then use linear 
algebra to recover the coefficients b n . 

Remark 3.3. In practice, it is better to work not with f(z) but instead the normalized 
function 

/norm(z)=/(z)(Im Z )*/ 2 

since then we have the transformation formula 

/normfe) = f(gz)(lmg Z f 2 = j(g,z) k f(z) (^^2)^ ' = (|J^T ) Vm-mW 

' " (12) 

and the automorphy factor has absolute value 1 . In particular, f notm is bounded on 
D (since it is bounded on the fundamental domain D), so working with / nor m yields 
better numerical stability if p is large, since then the contribution from the automor- 
phy factor may otherwise be quite large. Although /norm is no longer holomorphic, 
this does not affect the above method in any way. 



Using the Cauchy integral formula 



An alternative way to obtain the coefficients b„ is to use the Cauchy integral formula: 
for n > 0, we have 



bn = 7T— / T7T^ TT dw 

2niJ w n+l {\-w) k 
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with the integral a simple contour around 0. If we take this contour to be a circle of 
radius R> p, then again using automorphy, Cauchy's integral is equivalent to one 
evaluated along a path in D where the approximation f(z) w /n{z) holds; then using 
techniques of numerical integration we obtain a nontrivial linear relation among the 
coefficients b n . 

The simplest version of this technique arises by using simple Riemann summa- 
tion. Letting w = pe' e and dw = ipe' e d9, we have 

1 [*< fjpf) 

" 2%k {pe w ) n {\-pe w f ' 

again breaking up [0,2k] into Q e Z>3 intervals, letting w m = pe 27Cm, /2 an( j w ' m — 
g m w,„ with w' m e D, and z' m = z(p; w' m ) and z m = z(p; w m ), we obtain 

b _ I y f( Z >«) ~ 1 V J(8m,Z m y k f N (z' m ) 

and expanding /at(z) we obtain 

N 

b n ^Y, K nrbr (13) 

r=0 



where ("c" for Cauchy) 



K c — — V j(8m,z m ) k . , y 



with an error in the approximation that is small if Q is large. The matrix K c with 
entries with rows indexed by n = 0, . . . , and columns indexed by r = 0, . . . , N, 
is obtained from (13) by a matrix multiplication: 

QK C =JW' (14) 

where J is the matrix with entries 



j {gm j 2m 



^fim „ /, \ j, (15) 



w n m {l-w m ) 

with <n<N and 1 <m<Q and W 7 is the Vandermonde matrix with entries 

W' = (vJ Y 



with 1 < m < Q and < r <N. The matrices J and W' are fast to compute, so com- 
puting K c comes essentially at the cost of one matrix multiplication. The column 
vector b with entries b n satisfies K c b w b, so b is approximately in the kernel of 
K c -l. 
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Remark 3.4. For large values of n, the term w" n in the denominator of (15) dominates 
and creates numerical instability. Therefore, instead of solving for the coefficients 
b n we write b' n — b n p n and 

CO 

/(z) = (i-w)*i»/pr 

n=0 

and solve for the coefficients b' n . This replaces w m by w m /p with absolute value 
1. We then further scale the relation (13) by 1 /{K„ n — 1) so that the coefficient of 
b n in the relation is equal to 1; in practice, we observe that the largest entry in 
absolute value in the matrix of relations occurs along the diagonal, yielding quite 
good numerical stability. 

Better still is to use Simpson's rule: Break up the interval [0,2ft] into 2Q+ 1 e 
Z>3 intervals of equal length. With the same notation as above we then have 

, _ J_ y ( j{g2m-i,Z2m-iY k f N {w' 2m _ x ) ^ jfem^m)'* '/n^J 
" " 6G h V W 2 m -1 (1 - ™2,n-lT W" 2m (l - W2m f 

+ j(g2m+\,Z2m+\Y k fN(w2 m+ l) \ 
W n 2m+x {l-W2m+\) k )' 

(16) 

We obtain relations analogous to (13) and (14) in a similar way, the latter as a sum 
over three factored matrices. Then we have 

N 

b„~Y. L nr b r 

where 

TC 1 V 1 ( j(82m-l,Z2m-iyk , v . }(g2m,Z2mY k , 

j{g2m+\,Z2m+\)~k , \ 

Then, letting L c denote the matrix formed by L L nr with rows indexed by n and 
columns indexed by r, we have that Lfb w b, where b is the column vector hav- 
ing as its entries the b n , and hence b is in the numerical kernel of L c — 1. 

Simpson's rule is fast to evaluate and gives quite accurate results and so is quite 
suitable for our purposes. In very high precision, one could instead use more ad- 
vanced techniques for numerical integration; each integral in D can be broken up 
into a finite sum with contours given by geodesies. 

Remark 3.5. The coefficients b n fovn>N are approximately determined by the 
coefficients b n for n < N by the relation (13). In this way, they can be computed 
using integration and without any further linear algebra step. 
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Hecke operators 

In the special situation where F is a congruence arithmetic Fuchsian group, we can 
impose additional linear relations on eigenforms by using the action of the Hecke 
operators, as follows. 

We use the notation introduced in Section 1 . Suppose that / is an eigenform for 
F (still of weight k) with F of level 91. Let p be a nonzero prime of Zp with p \ £>9T. 

Using methods of Greenberg and the first author [21], we can compute for every 
prime p \ DVl the Hecke eigenvalue a p of an eigenform / using explicit methods 
in group cohomology. In particular, for each p we compute elements %\ , . . . , n q € 
PSL 2 (£) C PSL 2 (M) with q = Np + 1 so that the action of the Hecke operator T p is 
given by 

(Tpf)(z) = Y,j{m,z)- k f{mz) = a p f(z). (17) 

Taking z = p, for example, and writing w p j = w{%ip) for i = 1,.. . ,q and w' p i — 
giWi G D as before, we obtain 

« P /(o) ~ij(n l , P )- k j( 8l , Wp ,y k f(z P , l )- 

i=l 

Expanding / as a series in w as before, we find 

a„/(0)«£O» (18) 

n=0 

where ("h" for Hecke) 

K h n =im,pr k jOsi,*pj- k (i-^) k WJ n - 

i=l 

One could equally well consider the relations induced by plugging in other values 
of z, but in our experience are not especially numerically stable; rather, we use a few 
Hecke operators to isolate a one-dimensional subspace in combination with those 
relations coming from modularity. 

Remark 3.6. Alternatively, having computed the space Mfc(F), one could turn this 
idea around and use the above relations to compute the action of the Hecke operators 
purely analytically! For each / in a basis for M^T), we evaluate T p f at enough 
points to write T p f in terms of the basis, thereby giving the action of T p on M^(F). 
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Thus far, to encode the desired relationships on the coefficients, we have used the 
action of the group (i.e. automorphy) and Hecke operators. We may obtain further 
relationships obtained from the Shimura-Maass derivatives. The following lemma 
describes the action of d on power series expansions (see also Datskovsky and 
Guerzhoy [14, Proposition 2]). 

Lemma 3.7. Let f e M k (r) and p € !K. Then 
(d m f)(z) = T h (OT + ; )fc -^- r(w) (l-w)^ Y (dn+rf){p \ -Any)^ 



where y = Imp and 



Proof. The result holds for m — 1 by Lemma 2.2. Now suppose that the Lemma 
holds up to some positive integer m > 1 . We have 

(d m+1 f)(z) = d k+2m (d m f)(z) =(~- id m f)iz). 



Inidz 4^Im(z) 

Hence, after substituting for {d m f){z) and observing that dw/dz = (1 — w) 2 / (2iy), 
we obtain an expression for (<9 m+1 /)(z), in which we observe that for all n € Z>o 



1 ds n {w) = J_ y( _ 1)n - t (n\ (1-w)' 
2ni dz 4n^ Q [ ' \t J y>Im(z)"-> 



t=0 

Utilizing this in conjunction with 



f(l-vv) n-t 



Im(z) 



1 ds m - r (w) (k + 2r)s m - r (w)(l -w) (k + 2m)s m _ r (w) 



2ni dz 4-ny Anlm(z) 

_ (k + m + r)s m+ i- r (w) 
4n 

we obtain the result for d' n+l , from which the result follows by induction. 

For an eigenform / £ M^{r), we can further use the Hecke operators in conjunc- 
tion with the Shimura-Maass derivatives. 

Proposition 3.8 (Beyerl, James, Trentacose, Xue [7]). Let f e M k (r). Then d™(f) 
is a Hecke eigenform if and only iff is and eigenform; if so, and a„ is the eigenvalue 
ofT n associated to f, then the eigenvalue ofT n associated to d™(f) is n m a n . 
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Following the previous subsections, we assemble linear relations into a matrix A 
with M rows and N +1 columns such that Ab w 0. We now turn to compute the 
numerical kernel of A. We may assume that M > N but do not necessarily require 
that M = N. 

Suppose first we are in the case where the space of forms of interest is one- 
dimensional. This happens when, for example, X — r\K has genus g = 1 and k = 2; 
one can also arrange for this as in the previous section by adding linear relations 
coming from Hecke operators. Then since always belongs to the numerical ker- 
nel, we can dehomogenize by setting bo = 0; then letting A 1 be the M x N matrix 
consisting of all columns of A but the first column, let b' be the column vector with 
unknown entries b\,..., bn, and let a be the first column of A. Then A'b' — —a, and 
we can apply any numerically stable linear algebra method, an LU decomposition 
for example, to find a (least-squares) solution. 

There is a more general method to compute the entire numerical kernel of A 
which yields more information about the numerical stability of the computation. 
Namely, we compute the singular value decomposition (SVD) of the matrix A, writ- 
ing 

A = USV* 

where U and V are M x M and (N+ 1) x (N+ 1) unitary matrices and S is diagonal, 
where V* denotes the conjugate transpose of V. The diagonal entries of the matrix S 
are the singular values of A, the square roots of the eigenvalues of A* A, and may be 
taken to occur in decreasing magnitude. Singular values that are approximately zero 
correspond to column vectors of V that are in the numerical kernel of A, and one 
expects to have found a good quality numerical kernel if the other singular values 
are not too small. 



Confirming the output 

We have already mentioned several ways to confirm that the output looks correct. 
The first is to simply decrease e and see if the coefficients b n converge. The second 
is to look at the singular values to see that the approximately nonzero eigenvalues 
are sufficiently large (or better yet, that the dimension of the numerical kernel is 
equal to the dimension of the space M^r), when it can be computed using other 
formulas). 

More seriously, we can also check that the modularity relations (10) hold for a 
point w E D with \w\ < p but w D. (Such points always exist as the fundamental 
domain D is hyperbolically convex.) This test is quick and already convincingly 
shows that the computed expansion transforms like a modular form of weight k. 

Finally, when / is an eigenform for a congruence group T, we can check that 
/ is indeed numerically an eigenform (with the right eigenvalues) and that the co- 
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efficients, when normalized as (5), appear to be algebraic using the LLL-algorithm 
[26]. 



4 Results 

In this section, we present four examples to demonstrate our method. 



Example 1 

We begin by computing with a classical modular form so that we can use its q- 
expansion to verify that the power series expansion is correct. However, our method 
does not work well in this non-cocompact situation, so we prepare ourselves for 
poor accuracy. 

Let / e £2(^0(1 1)) be the unique normalized eigenform of weight 2 and level 11, 
defined by 



where q = e 2mz . We choose a CM point (Heegner point) for lo(ll) for the field 
K = Q(V^7) with absolute discriminant d = 7, namely p = (-9 + \/ = 7) /22. We 
find the power series expansion (4) written (5) directly using Lemmas 2.2 and 2.4: 
we obtain 



f(z) = qW <f?( l <7 11 ") 2 = 1 2 <? 2 - <? 3 + 2 ? 4 + - = E a n4 



n 




m = a - w) 2 £ b„w- = f(p)(i wf £ %©w) 



n=0 n=0 




where 




and 
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Fig. 1: A fundamental domain for the modular curve Xq(1 1) 



, (d-\ \ '/ 2 'M 

l / r(i/7)r(2/7)r(4/7) \ 1/2 

" \T (3 /7)r (5 /7)r (6/7) J 



= 0.5004912..., 



and — a/3 + 4^7—7 = —2.6457513 . . . +2z. The coefficients c„ are found always to 
be integers, and so are obtained by rounding the result obtained by numerical ap- 
proximation; however, we do not know if this integrality result is a theorem, and 
consequently this expansion for / is for the moment only an experimental observa- 
tion. 

We now apply our method and compare. We first compute a fundamental domain 
for Jo (11) with center p, shown in Figure 1. 

Now we take e = 10~ 20 . In this situation, the fundamental domain is not compact, 
so we must adapt our method. We choose a radius p for which the image under 
io(ll) to the fundamental domain will yield nontrivial relations; we choose p = 
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0.85. Then by (8), we may take N = 300 (adding some for good measure) for the 
estimate f(z) s=s /n(z) whenever \w(z)\ < p. 

We compute the relations (13) from the Cauchy integral formula for radius R, and 
we add the relations (18) coming from Hecke operators for the primes 2,3,5,7, 13. 
(The fact that the image of a point in the fundamental domain under the Hecke 
operators can escape to a point in the fundamental domain but with radius > R 
precludes the further use of higher Hecke operators.) 

We compute the SVD of this matrix of relations and find the smallest singular 
value to be < e; the next largest singular value is only 4 times as large, so also 
arguably negligible. In any case, the associated element of the numerical kernel 
corresponding to the smallest singular value yields a series (1 — w) k Y%=ob n w n with 
R n \b„ — b n \ < 10~ 9 for the first 10 coefficients but increasingly inaccurate for n 
moderate to large. Given that we expected unsatisfactory numerical results, this test 
at least indicates that our matrix of relations (coming from both the Cauchy integral 
formula and Hecke operators) passes one sanity check. 



Example 2 

For a comparison, we now consider the well-studied example arising from the 
(2,4, 6) -triangle group associated to the quaternion algebra of discriminant 6, refer- 
enced in the introduction. 

We follow Bayer [2], Bayer and Travesa [3, 4], and Baba and Granath [1]; for 
further detail, we refer to these articles. 
/3,-l\ 

Let B = — be the quaternion algebra of discriminant 6 over the rationals 



Q, so that a 2 = 3, j3 2 = -1, and pa = -a/3. A maximal order C B is given by 

= Z®aZ®j3Z©5Z 
where 8 = (1 + a + j8 + aj8)/2. We have the splitting 



loo :B^M 2 { 



fV3 WO 1 



a > P "{0 -V3, 



-1 



Let r = loo(Oj X )/{±l}. Then T is a group of signature (0;2,2,3,3) which is an 
index 4 normal subgroup in the triangle group 

4(2,4,6) = (72,74,76 I 72 = H = l£ = 727476 = 1). 

The point 

V6-V2 

P= o 1 
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Fig. 2: A fundamental domain for the Shimura curve X associated to a maximal 
order in the quaternion algebra of discriminant 6 over Q 



is a CM point of absolute discriminant d = 24. The Dirichlet fundamental domain 
D{p) is then given in Figure 2. (A slightly different fundamental domain, empha- 
sizing the relationship to the triangle group 4(2,4,6), is considered by the other 
authors.) We have p = 0.447213 .... 

We consider the space S^r) of modular forms for F of weight 4; by Riemann- 
Roch, this space has dimension 1 and so is generated by a modular form P e Sn{r) 
with expansion 

P(z) = (V=3(2 V3)Q 4 ) (1 - w) 4 (1 + A i. ( 0W )2 _ « J_ (0w) 4 + 

555 1 _ . fi 57165 1 ,^ .» \ 

W @W) + ^8! (0W) + -J 

(20) 

where @ = -4^X2 2 and 
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, (d-x \ '/ 2 »M 

Q = ^= r\r{j/d) { - d/j) =0.321211772390... 
V2dn \j = \ J 

with d = 24. 

Remark 4. 1. The algebraic factor y = lm(p) = \J2{ y/3 — 1 ) /2 need not be multipled 
in the period 0, as expected by (5), to obtain coefficients defined over Q. 

We present the algorithm using the formulas obtained from Simpson's rule. Let 
K c be the (N + 1) x (N + 1) matrix obtained from Simpson's rule. We then take 
a square N x N submatrix of full rank and set Icq to be the removed row vector to 
obtain a system Ab = ko; we wish to solve for b. An LU decomposition and back 
substitution is then used to obtain coefficients b n normalized so that bo = 1. Let 
£>«,exact denote the exact coefficients in (20) above, renormalized so that £>o,exact = 1 ■ 
We set Q = 2N as this choice provides a good approximation to the integral. Then 
we obtain the following errors. 



Table 1 : Example 2 Results 



N 


P|&l,exact-&l| 


max„p" fc„,exact-&n 


35 


io- 13 


io- 13 


70 


io- 23 


10~ 22 


140 


io- 47 


io- 47 



We find that with a fixed radius, we obtain better precision in the answer as Af 
is increased. That is, the results are completely determined, for a fixed radius and 
precision, by the degree of the approximating polynomial. 

Remark 4.2. In contrast to the non-cocompact case, there is no gain in the accuracy 
of the results from taking a larger radius. Hence it suffices to fix the radius p. 



Example 3 

Next, we work with an arithmetic group over a totally real field, and compute the 

image of a CM point under the Shimura curve parametrization of an elliptic curve. 

Let F = Q(a) = Q(v / 5) where a 2 + a — 1 =0, and let Zp be its ring of integers. Let 

p = (5a + 2), so Np =31. Let B be the quaternion algebra ramified at p and the real 

r /a,5fl + 2\ 
place sending v 5 to its positive real root: we take B = I — I . We consider 

F E embedded by the second real place, so a = (1 -y/5)/2 = -0.618033.... 
A maximal order C B is given by 
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,o ™ n, a + aa-p „ (a-l)+aa-aB 

= Z F ®Z F a®Z F — -®Z F ^ 

Let loo be the splitting at the other real place given by 

„ / Ja\ f^5a + 2 



y ' V -V5a + 2 / 

LetT = ^(Of )/{±l} C PSL 2 (M). Then T has hyperbolic area 1 (normalized so 
that an ideal triangle has area 1 /2) and signature (1;2 2 ), so X = F\!K can be given 
the structure of a compact Riemann surface of genus 1. Consequently, the space 
52(F) of modular forms on F of weight 2 is 1 -dimensional, and it is this space that 
we will compute. 

The field K = F(\/=7) embeds in with 

1 5a+10 a + 2 p 3a-5 . m 
M = "2 2~«-— ^ + — 

satisfying fx 2 + fi +2 = and Z F [jit] = Z# the maximal order with class number 1. 
We take p = -3.1653 . . . + 1.41783 . . . ^/^T e % to be the fixed point of ju, a CM 
point of discriminant —7. 

We compute a Dirichlet fundamental domain D(p) for F as in Figure 3. 

We have p = 0.71807 ... so for e = 10~ 20 we take N = 150, and we take R = 
0.8. As 52(F) is 1 -dimensional, we use only the relations coming from the Cauchy 
integral formula (13) and reserve the relations from the Hecke operators as a check. 
The (N+ 1) x (N+ l)-matrix has largest singular value 4.01413 . . . and one singular 
value which is < e — the next largest singular value is 0.499377 . . ., showing that the 
numerical kernel is one-dimensional. 

Computing this kernel, we find numerically that 

2/1 { c\ \ 70a + 114 2 8064a + 13038 3 



(21) 



f(z) = (l-w) z [l + (0w) {0wY - (0w) 

174888a + 282972,^ n4 13266960a + 21466440 ln n5 
4j (®vv) 4 {QwY- 

1826784288a + 2955799224 _ n6 

gj N- 

2388004416a + 3863871648 _ x7 

(@w) + . 



7! 
where 

= 0.046218579529208499918 ... - 0.075987317531832568351 . 



is a period akin to (19) and related to the CM abelian variety associated to the 
point p. (We do not have a good way to scale the function /, so we simply choose 
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Fig. 3: A fundamental domain for the Shimura curve X associated to a maximal 
order in the quaternion algebra ramified at a prime above 3 1 and the first real place 
overF =Q(\/5) 

f(p) = 1 .) The apparent integrality of the numerators allows the LLL-algorithm to 
identify these algebraic numbers even with such low precision. 

We then also numerically verify (to precision e) that this function is an eigen- 
function for the Hecke operators for all primes [ with M < 30. This gives convincing 
evidence that our function is correct. 

We further compute the other embedding of this form by repeating the above 
computation with conjugate data: we take an algebra ramified at p and the other real 
place. We find the following fundamental domain in Figure 4. 

A similar computation yields a matrix of relations with a one-dimensional nu- 
merical kernel whose corresponding expansion has coefficients which agree with 
the conjugates of those in (21) under a i-> — (a + 1), the generator of the group 
Gal(Q(x/5)/Q). 
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Fig. 4: A fundamental domain for the Shimura curve X associated to a maximal 
order in the quaternion algebra ramified at a prime above 3 1 and the second real 
place over F = Q(\/5) 



Next, we identify the equation of the Jacobian J of the curve X by computing 
the associated periods. We first identify the group F using the sidepairing relations 
coming from the computation of D(p) [47]: 

r = <7i, 72, 8 l ,8z\$i= Si = yf 1 y^ 1 5, y, ji&i = 1) 



where 



a + 2 2a + 3 a+l „ 
7i = ^ ^a + — a/3 

2a + 3 7a + 10 a + 2 . f . „ 

y 2 = — — + — - — a + — — j3-(3a + 5)aj3 



generate the free part of the maximal abelian quotient of F. The elements y\,f2 
identify vertices: V2 i-> V5 = 71 (V2) and vs V2 = 72(vg). Therefore, we compute 
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two independent periods 0)i , a>2 (up to scaling) 

n _ dw ( N 



Jv 2 (1-w) 2 \„ i To"+ 1 



"5 

= -0.654017... +0.397799... i 

"2 



dw 



0)2= / /(z W - ., =0.952307... +0.829145 
We then compute the j-invariant 

1 188961 1722383394a + 86293850621 19691 



j{0h/(Sh) = -18733.423... = -- 



31 8 



and by computing the values of the Eisenstein series £4 (©1/0)2) and E^Obi/obi), or 
by twisting a curve with j-invariant 7(0)1/0)2), we identify the elliptic curve J as 

y 2 +xy-ay = x i - (a - l)x 2 - (31a + 75)x - (141a + 303). 

We did this by recognizing this algebraic number using LLL and then confirming by 
computing the conjugate of j under Gal(Q(v / 5)/Q) as above and recognizing the 
trace and norm as rational numbers using simple continued fractions. Happily, we 
see that this curve also appears in the tables of elliptic curves over Q(\/5) computed 
[15,42]. 

Finally, we compute the image on J of a degree zero divisor on X. The fixed 
points w\ , W2 of the two elliptic generators 81 and ^2 are CM points of discriminant 
—4. Let K = F(i) and consider the image of [wi] — [w2\ on J given by the Abel- 
Jacobi map as 

r w 2 dw 

/ /(z)t; v, = -0.177051... -0.291088... i (mod A) 

Jwj (l — wY 

where A = Zo)i + Zo)2 is the period lattice of J. Evaluating the elliptic exponential, 
we find the point 

(-10.503797 . . . ,5.560915 ... - 44. 133005 . . . i) e J{C) 

which matches to the precision computed e = 10~ 20 the point 

-81 fl - 118 (358 fl+ 1191)/ + (194 fl + 236)\ e z 
16 64 / 

which generates J(K)/J(K) t0!:s . On the other hand, the Mordell-Weil group of J 
over F is in fact torsion, with J(F) = Z/2Z, generated by the point ((28a + 
27)/4,-(24a + 27)/8). 
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Example 4 



To conclude, we compute the equation of a genus 2 Shimura curve over a totally real 
field. (There is nothing special about this curve; we chose it because only it seemed 
to provide a nice example.) 



Let F 



(a) 



13) where w 



w - 



3 = 0, and let Z F = Z[w] be its 



ring of integers. Let p 

, ,'w-2,2w-l 
bra B = 



(2w — 1) = (yl3). We consider the quaternion alge- 
which is ramified at p and the real place sending a 

r i 

(1 - \/l3)/2 = - 1 .302775 .... We consider F R embedded by the first real place, 
so w = (l+V\3)/2 = 2.302775... . 
A maximal order C B is given by 



= Z/r®Z F a©Z F 



wa + j3 (w + l) + a/3 



Let loo be the splitting 



loo:B^M 2 ( 



a,j3 !->• 













-\/2v 



Let r = loo(0f ) /{±1} C PSL 2 (M). Then T has hyperbolic area 2 and signature 
(2;—), so X = r\K can be given the structure of a compact Riemann surface of 
genus 2 and 5 2 (F) is 2-dimensional. We compute [21] that this space is irreducible 
as a Hecke module, represented by a constituent eigenform / with the following 
eigenvalues: 



P 

Np 



«*(/) 



(w) (w-1) (2) (2w-l) (w + 4) (w-5) (3w + l) (3w-4) 
3 3 4 13 17 17 23 23 



-e+1 



-2 -1 2e-4 -2e-2 -2e + 4 2e + 2 



Here, the element e E Q satisfies e 2 — e — 5 = and the Hecke eigenvalue field 
£ = Q(e) has discriminant 21. Let x £ Gal(£/Q) represent the nontrivial element; 
then 5 2 (-T) is spanned by / and t(/), where a p (t(/)) = t(o p (/)). 

We compute a Dirichlet fundamental domain D(p) for F as in Figure 5. 

We have p = 0.88561 1 ... so for e = lO" 20 we take N = 600. For the form / and 
its conjugate t(/), we add the implied relations for the Hecke operators to those 
from the Cauchy integral. In each case, we find one singular value < e and the next 
largest singular value to be > 0.2, so we have indeed found two one-dimensional 
kernels. 

We first consider g = f + t(/), and find numerically that 



Computing power series expansions of modular forms 



27 




Fig. 5: A fundamental domain for the Shimura curve X associated to a maximal 
order in the quaternion algebra ramified at a prime above 13 and the second real 
place over F = Q(vT3) 



, , „ , n . 28496,^ , 2 14796288 _ 
g(z) = (1 - wf ( 1 + 117(0w) + — (0w) 2 (0w) 3 

3287025664,^ , 4 1142501867520,^ n5 
4j (0w) (0w) 5 

116349452943360, xfi 420556103693107200, x7 
gj (0w) (0 w) 1 

where 

= 0.008587333922373292375569160851 .... 



(22) 
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We see in particular that the coefficients are rational numbers. 
We then compute h= (/ — t(/))/\/21 and find 




i&w) 1 + ... 



) 



(23) 



Then, from this basis of differentials, we can use the standard Riemann-Roch 
basis argument to compute an equation for the curve X (see e.g. Galbraith [20, §4]). 
Let x = g/h and y — x' /h, where ' denotes the derivative with respect to ©w; then 
since h has a zero of order 1 at p, the function x has a pole of order 1 , and similarly 
y has a pole of order 3, therefore by Riemann-Roch we obtain an equation y 2 = q{x) 
with^(x) e Q[x]: 



and this reduced model has discriminant 2995508600908518877 = 7 10 13 9 . 

The curve X over F, however, has good reduction away from £> = (Vl3) [10]. 
Indeed, since we don't know the precise period to normalize the forms, we may 
inadvertently introduce a factor coming from the field K = F(y/—7) of CM of the 
point p at which we have expanded these functions as a power series. In this case, 
our model is off by a quadratic twist coming from the CM field K = F{\/^1): 
twisting we obtain 



The only automorphism of X is the hyperelliptic involution, so Aut(X) = Z/2Z. 

Since F is Galois over Q, has narrow class number 1, and the discriminant S of 
B is invariant under Gal(F/Q), it follows from the analysis of Doi and Naganuma 
[17] that the field of moduli of the curve X is Q. This does not, however, imply that 
X admits a model over Q: there is an obstruction in general to a curve of genus 2 to 
having a model over its field of moduli (Mestre's obstruction [27], see also Cardona 
and Quer [11]). This obstruction apparently vanishes forX; it would be interesting 
to understand this phenomenon more generally. 

Finally, this method really produces numerically the canonical model of X (equa- 
tion (24) considered over the reflex field F), in the sense of Shimura [35]. In this 
way, we have answered "intriguing question" of Elkies [18, 5.5]. 



y 2 = x 6 + 1950x 5 + 828919x 4 - 122128188x 3 + 3024544159* 2 
-29677133122*+ 107045514121. 



This curve has a reduced model 



y 2 + y = lx 5 - 84x 4 + 1 19x 3 + 749x 2 + 938x + 390 



X : y 2 +y = x 5 + I2x 4 + 17x 3 - 107x 2 + 134x - 56. 



(24) 
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